home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / nurb_polyg / NurbSubdiv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  17.4 KB  |  643 lines
  1. /*
  2.  * NurbSubdiv.c - Perform adaptive subdivision on a NURB surface.
  3.  *
  4.  * John Peterson
  5.  *
  6.  */
  7.  
  8. #include <stdlib.h>
  9. #include <stdio.h>
  10. #include <math.h>
  11.  
  12. #include "nurbs.h"
  13. #include "drawing.h"
  14.  
  15. #define EPSILON 0.0000001   /* Used to determine when things are too small. */
  16. #define DIVPT( p, dn ) { ((p).x) /= (dn); ((p).y) /= (dn); ((p).z) /= (dn); }
  17.  
  18. double SubdivTolerance;        /* Size, in pixels of facets produced */
  19.  
  20. #define maxV(surf) ((surf)->numV-1L)
  21. #define maxU(surf) ((surf)->numU-1L)
  22.  
  23. /*
  24.  * Split a knot vector at the center, by adding multiplicity k knots near
  25.  * the middle of the parameter range.  Tries to start with an existing knot,
  26.  * but will add a new knot value if there's nothing in "the middle" (e.g.,
  27.  * a Bezier curve).
  28.  */
  29. static long
  30. SplitKV( double * srckv,
  31.      double ** destkv,
  32.      long * splitPt,    /* Where the knot interval is split */
  33.      long m, long k )
  34. {
  35.     long i, last;
  36.     long middex, extra, same;    /* "middex" ==> "middle index" */
  37.     double midVal;
  38.  
  39.     extra = 0L;
  40.     last = (m + k);
  41.  
  42.     middex = last / 2;
  43.     midVal = srckv[middex];
  44.  
  45.     /* Search forward and backward to see if multiple knot is already there */
  46.  
  47.     i = middex+1L;
  48.     same = 1L;
  49.     while ((i < last) && (srckv[i] == midVal)) {
  50.     i++;
  51.     same++;
  52.     }
  53.  
  54.     i = middex-1L;
  55.     while ((i > 0L) && (srckv[i] == midVal)) {
  56.     i--;
  57.     middex--;        /* middex is start of multiple knot */
  58.     same++;
  59.     }
  60.  
  61.     if (i <= 0L)        /* No knot in middle, must create it */
  62.     {
  63.     midVal = (srckv[0L] + srckv[last]) / 2.0;
  64.     middex = last / 2L;
  65.     while (srckv[middex + 1L] < midVal)
  66.         middex++;
  67.     same = 0L;
  68.     }
  69.  
  70.     extra = k - same;
  71.     CHECK( *destkv = (double *) malloc( (long) (sizeof( double )
  72.                     * (m+k+extra+1L) ) ) );
  73.  
  74.     if (same < k)        /* Must add knots */
  75.     {
  76.     for (i = 0L; i <= middex; i++)
  77.         (*destkv)[i] = srckv[i];
  78.  
  79.     for (i = middex+1L; i <= middex+extra; i++)
  80.         (*destkv)[i] = midVal;
  81.  
  82.     for (i = middex + k - same + 1L; i <= m + k + extra; i++)
  83.         (*destkv)[i] = srckv[i - extra];
  84.     }
  85.     else
  86.     {
  87.     for (i = 0L; i <= m + k; i++)
  88.     (*destkv)[i] = srckv[i];
  89.     }
  90.  
  91.     *splitPt = (extra < k) ? middex - 1L : middex;
  92.     return( extra );
  93. }
  94.  
  95. /*
  96.  * Given a line defined by firstPt and lastPt, project midPt onto
  97.  * that line.  Used for fixing "cracks".
  98.  */
  99. static void
  100. ProjectToLine( Point3 * firstPt, Point3 * lastPt, Point3 * midPt )
  101. {
  102.     Point3 base, v0, vm;
  103.     double fraction, denom;
  104.  
  105.     base = *firstPt;
  106.  
  107.     (void) V3Sub( lastPt, &base, &v0 );
  108.     (void) V3Sub( midPt, &base, &vm );
  109.  
  110.     denom = V3SquaredLength( &v0 );
  111.     fraction = (denom == 0.0) ? 0.0 : (V3Dot( &v0, &vm ) / denom);
  112.  
  113.     midPt->x = base.x + fraction * v0.x;
  114.     midPt->y = base.y + fraction * v0.y;
  115.     midPt->z = base.z + fraction * v0.z;
  116. }
  117.  
  118. /*
  119.  * If a normal has collapsed to zero (normLen == 0.0) then try
  120.  * and fix it by looking at its neighbors.  If all the neighbors
  121.  * are sick, then re-compute them from the plane they form.
  122.  * If that fails too, then we give up...
  123.  */
  124. static void
  125. FixNormals( SurfSample * s0, SurfSample * s1, SurfSample * s2 )
  126. {
  127.     Boolean goodnorm;
  128.     long i, j, ok;
  129.     double dist;
  130.     SurfSample * V[3];
  131.     Point3 norm;
  132.  
  133.     V[0] = s0; V[1] = s1; V[2] = s2;
  134.  
  135.     /* Find a reasonable normal */
  136.     for (ok = 0, goodnorm = FALSE;
  137.     (ok < 3L) && !(goodnorm = (V[ok]->normLen > 0.0)); ok++);
  138.  
  139.     if (! goodnorm)    /* All provided normals are zilch, try and invent one */
  140.     {
  141.     norm.x = 0.0; norm.y = 0.0; norm.z = 0.0;
  142.  
  143.     for (i = 0; i < 3L; i++)
  144.     {
  145.         j = (i + 1L) % 3L;
  146.         norm.x += (V[i]->point.y - V[j]->point.y) * (V[i]->point.z + V[j]->point.z);
  147.         norm.y += (V[i]->point.z - V[j]->point.z) * (V[i]->point.x + V[j]->point.x);
  148.         norm.z += (V[i]->point.x - V[j]->point.x) * (V[i]->point.y + V[j]->point.y);
  149.     }
  150.     dist = V3Length( &norm );
  151.     if (dist == 0.0)
  152.         return;        /* This sucker's hopeless... */
  153.  
  154.     DIVPT( norm, dist );
  155.  
  156.     for (i = 0; i < 3; i++)
  157.     {
  158.         V[i]->normal = norm;
  159.         V[i]->normLen = dist;
  160.     }
  161.     }
  162.     else        /* Replace a sick normal with a healthy one nearby */
  163.     {
  164.     for (i = 0; i < 3; i++)
  165.         if ((i != ok) && (V[i]->normLen == 0.0))
  166.         V[i]->normal = V[ok]->normal;
  167.     }
  168.     return;
  169. }
  170.  
  171. /*
  172.  * Normalize the normal in a sample.  If it's degenerate,
  173.  * flag it as such by setting the normLen to 0.0
  174.  */
  175. static void
  176. AdjustNormal( SurfSample * samp )
  177. {
  178.     /* If it's not degenerate, do the normalization now */
  179.     samp->normLen = V3Length( &(samp->normal) );
  180.  
  181.     if (samp->normLen < EPSILON)
  182.     samp->normLen = 0.0;
  183.     else
  184.     DIVPT( (samp->normal), samp->normLen );
  185. }
  186.  
  187. /*
  188.  * Compute the normal of a corner point of a mesh.  The
  189.  * base is the value of the point at the corner, indU and indV
  190.  * are the mesh indices of that point (either 0 or numU|numV).
  191.  */
  192. static void
  193. GetNormal( NurbSurface * n, long indV, long indU )
  194. {
  195.     Point3 tmpL, tmpR;    /* "Left" and "Right" of the base point */
  196.     SurfSample * crnr;
  197.  
  198.     if ( (indU && indV) || ((! indU) && (!indV)) )
  199.     {
  200.     if (indU)
  201.         crnr = &(n->cnn);
  202.     else
  203.         crnr = &(n->c00);
  204.     DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), &tmpL );
  205.     DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), &tmpR );
  206.     }
  207.     else
  208.     {
  209.     if (indU)
  210.         crnr = &(n->c0n);
  211.     else
  212.         crnr = &(n->cn0);
  213.     DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), &tmpR );
  214.     DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), &tmpL );
  215.     }
  216.  
  217.     (void) V3Sub( &tmpL, &(crnr->point), &tmpL );
  218.     (void) V3Sub( &tmpR, &(crnr->point), &tmpR );
  219.     (void) V3Cross( &tmpL, &tmpR, &(crnr->normal) );
  220.     AdjustNormal( crnr );
  221. }
  222.  
  223. /*
  224.  * Build the new corners in the two new surfaces, computing both
  225.  * point on the surface along with the normal.    Prevent cracks that may occur.
  226.  */
  227. static void
  228. MakeNewCorners( NurbSurface * parent,
  229.         NurbSurface * kid0,
  230.         NurbSurface * kid1,
  231.         Boolean dirflag )
  232. {
  233.     DIVW( &(kid0->points[maxV(kid0)][maxU(kid0)]), &(kid0->cnn.point) );
  234.     GetNormal( kid0, maxV(kid0), maxU(kid0) );
  235.  
  236.     if (dirflag)
  237.     {
  238.     kid0->strUn = FALSE;    /* Must re-test new edge straightness */
  239.  
  240.     DIVW( &(kid0->points[0L][maxU(kid0)]), &(kid0->c0n.point) );
  241.     GetNormal( kid0, 0L, maxU(kid0) );
  242.     /*
  243.      * Normals must be re-calculated for kid1 in case the surface
  244.      * was split at a c1 (or c0!) discontinutiy
  245.      */
  246.     kid1->c00.point = kid0->c0n.point;
  247.     GetNormal( kid1, 0L, 0L );
  248.     kid1->cn0.point = kid0->cnn.point;
  249.     GetNormal( kid1, maxV(kid1), 0L );
  250.  
  251.     /*
  252.      * Prevent cracks from forming by forcing the points on the seam to
  253.      * lie along any straight edges.  (Must do this BEFORE finding normals)
  254.      */
  255.     if (parent->strV0)
  256.         ProjectToLine( &(parent->c00.point),
  257.                &(parent->c0n.point),
  258.                &(kid0->c0n.point) );
  259.     if (parent->strVn)
  260.         ProjectToLine( &(parent->cn0.point),
  261.                &(parent->cnn.point),
  262.                &(kid0->cnn.point) );
  263.  
  264.     kid1->c00.point = kid0->c0n.point;
  265.     kid1->cn0.point = kid0->cnn.point;
  266.     kid1->strU0 = FALSE;
  267.     }
  268.     else
  269.     {
  270.     kid0->strVn = FALSE;
  271.  
  272.     DIVW( &(kid0->points[maxV(kid0)][0]), &(kid0->cn0.point) );
  273.     GetNormal( kid0, maxV(kid0), 0L );
  274.     kid1->c00.point = kid0->cn0.point;
  275.     GetNormal( kid1, 0L, 0L );
  276.     kid1->c0n.point = kid0->cnn.point;
  277.     GetNormal( kid1, 0L, maxU(kid1) );
  278.  
  279.     if (parent->strU0)
  280.         ProjectToLine( &(parent->c00.point),
  281.                &(parent->cn0.point),
  282.                &(kid0->cn0.point) );
  283.     if (parent->strUn)
  284.         ProjectToLine( &(parent->c0n.point),
  285.                &(parent->cnn.point),
  286.                &(kid0->cnn.point) );
  287.  
  288.     kid1->c00.point = kid0->cn0.point;
  289.     kid1->c0n.point = kid0->cnn.point;
  290.     kid1->strV0 = FALSE;
  291.     }
  292. }
  293.  
  294. /*
  295.  * Split a surface into two halves.  First inserts multiplicity k knots
  296.  * in the center of the parametric range.  After refinement, the two
  297.  * resulting surfaces are copied into separate data structures.     If the
  298.  * parent surface had straight edges, the points of the children are
  299.  * projected onto those edges.
  300.  */
  301. static void
  302. SplitSurface( NurbSurface * parent,
  303.           NurbSurface * kid0, NurbSurface * kid1,
  304.           Boolean dirflag )        /* If true subdivided in U, else in V */
  305. {
  306.     NurbSurface tmp;
  307.     double * newkv;
  308.     long i, j, splitPt;
  309.  
  310.     /*
  311.      * Add a multiplicty k knot to the knot vector in the direction
  312.      * specified by dirflag, and refine the surface.  This creates two
  313.      * adjacent surfaces with c0 discontinuity at the seam.
  314.      */
  315.  
  316.     tmp = *parent;    /* Copy order, # of points, etc. */
  317.     if (dirflag)
  318.     {
  319.     tmp.numU = parent->numU + SplitKV( parent->kvU,
  320.                         &newkv,
  321.                         &splitPt,
  322.                         maxU(parent),
  323.                         parent->orderU );
  324.     AllocNurb( &tmp, newkv, NULL );
  325.     for (i = 0L; i < tmp.numV + tmp.orderV; i++)
  326.         tmp.kvV[i] = parent->kvV[i];
  327.     }
  328.     else
  329.     {
  330.     tmp.numV = parent->numV + SplitKV( parent->kvV,
  331.                         &newkv,
  332.                         &splitPt,
  333.                         maxV(parent),
  334.                         parent->orderV );
  335.     AllocNurb( &tmp, NULL, newkv );
  336.     for (i = 0L; i < tmp.numU + tmp.orderU; i++)
  337.         tmp.kvU[i] = parent->kvU[i];
  338.     }
  339.     RefineSurface( parent, &tmp, dirflag );
  340.  
  341.     /*
  342.      * Build the two child surfaces, and copy the data from the refined
  343.      * version of the parent (tmp) into the two children
  344.      */
  345.  
  346.     /* First half */
  347.  
  348.     *kid0 = *parent;    /* copy various edge flags and orders */
  349.  
  350.     kid0->numU = dirflag ? splitPt+1L : parent->numU;
  351.     kid0->numV = dirflag ? parent->numV : splitPt+1L;
  352.     kid0->kvU = kid0->kvV = NULL;
  353.     kid0->points = NULL;
  354.     AllocNurb( kid0, NULL, NULL );
  355.  
  356.     for (i = 0L; i < kid0->numV; i++)    /* Copy the point and kv data */
  357.     for (j = 0L; j < kid0->numU; j++)
  358.         kid0->points[i][j] = tmp.points[i][j];
  359.     for (i = 0L; i < kid0->orderU + kid0->numU; i++)
  360.     kid0->kvU[i] = tmp.kvU[i];
  361.     for (i = 0L; i < kid0->orderV + kid0->numV; i++)
  362.     kid0->kvV[i] = tmp.kvV[i];
  363.  
  364.     /* Second half */
  365.  
  366.     splitPt++;
  367.     *kid1 = *parent;
  368.  
  369.     kid1->numU = dirflag ? tmp.numU - splitPt : parent->numU;
  370.     kid1->numV = dirflag ? parent->numV : tmp.numV - splitPt;
  371.     kid1->kvU = kid1->kvV = NULL;
  372.     kid1->points = NULL;
  373.     AllocNurb( kid1, NULL, NULL );
  374.  
  375.     for (i = 0L; i < kid1->numV; i++)    /* Copy the point and kv data */
  376.     for (j = 0L; j < kid1->numU; j++)
  377.         kid1->points[i][j]
  378.         = tmp.points[dirflag ? i: (i + splitPt) ][dirflag ? (j + splitPt) : j];
  379.     for (i = 0L; i < kid1->orderU + kid1->numU; i++)
  380.     kid1->kvU[i] = tmp.kvU[dirflag ? (i + splitPt) : i];
  381.     for (i = 0L; i < kid1->orderV + kid1->numV; i++)
  382.     kid1->kvV[i] = tmp.kvV[dirflag ? i : (i + splitPt)];
  383.  
  384.     /* Construct new corners on the boundry between the two kids */
  385.     MakeNewCorners( parent, kid0, kid1, dirflag );
  386.  
  387.     FreeNurb( &tmp );        /* Get rid of refined parent */
  388. }
  389.  
  390. /*
  391.  * Test if a particular row or column of control points in a mesh
  392.  * is "straight" with respect to a particular tolerance.  Returns true
  393.  * if it is.
  394.  */
  395.  
  396. #define GETPT( i )  (( dirflag ? &(n->points[crvInd][i]) : &(n->points[i][crvInd]) ))
  397.  
  398. static Boolean
  399. IsCurveStraight( NurbSurface * n,
  400.          double tolerance,
  401.          long crvInd,
  402.          Boolean dirflag )  /* If true, test in U direction, else test in V */
  403. {
  404.     Point3 p, vec, prod;
  405.     Point3 cp, e0;
  406.     long i, last;
  407.     double linelen, dist;
  408.  
  409.     /* Special case: lines are automatically straight. */
  410.     if ((dirflag ? n->numU : n->numV) == 2L)
  411.     return( TRUE );
  412.  
  413.     last = (dirflag ? n->numU : n->numV) - 1L;
  414.     ScreenProject( GETPT( 0L ), &e0 );
  415.  
  416.     /* Form an initial line to test the other points against (skiping degen lines) */
  417.  
  418.     linelen = 0.0;
  419.     for (i = last; (i > 0L) && (linelen < EPSILON); i--)
  420.     {
  421.     ScreenProject( GETPT( i ), &cp );
  422.     (void) V3Sub( &cp, &e0, &vec );
  423.  
  424.     linelen = sqrt( V3SquaredLength( &vec ) );
  425.     }
  426.  
  427.     DIVPT( vec, linelen );
  428.  
  429.     if (linelen > EPSILON)    /* If no non-degenerate lines found, it's all degen */
  430.     for (i = 1L; i <= last; i++)
  431.     {
  432.         /* The cross product of the vector defining the
  433.          * initial line with the vector of the current point
  434.          * gives the distance to the line. */
  435.         ScreenProject( GETPT( i ), &cp );
  436.         (void) V3Sub( &cp,&e0,&p );
  437.  
  438.         (void) V3Cross( &p, &vec, &prod );
  439.         dist = V3Length( &prod );
  440.  
  441.         if (dist > tolerance)
  442.         return( FALSE );
  443.     }
  444.  
  445.     return( TRUE );
  446. }
  447.  
  448. /*
  449.  * Check to see if a surface is flat.  Tests are only performed on edges and
  450.  * directions that aren't already straight.  If an edge is flagged as straight
  451.  * (from the parent surface) it is assumed it will stay that way.
  452.  */
  453. static Boolean
  454. TestFlat( NurbSurface * n, double tolerance )
  455. {
  456.     long i;
  457.     Boolean straight;
  458.     Point3 cp00, cp0n, cpn0, cpnn, planeEqn;
  459.     double dist,d ;
  460.  
  461.     /* Check edge straightness */
  462.  
  463.     if (! n->strU0)
  464.     n->strU0 = IsCurveStraight( n, tolerance, 0L, FALSE );
  465.     if (! n->strUn)
  466.     n->strUn = IsCurveStraight( n, tolerance, maxU(n), FALSE );
  467.     if (! n->strV0)
  468.     n->strV0 = IsCurveStraight( n, tolerance, 0L, TRUE );
  469.     if (! n->strVn)
  470.     n->strVn = IsCurveStraight( n, tolerance, maxV(n), TRUE );
  471.  
  472.     /* Test to make sure control points are straight in U and V */
  473.  
  474.     straight = TRUE;
  475.     if ( (! n->flatU) && (n->strV0) && (n->strVn) )
  476.     for (i = 1L;
  477.          (i < maxV(n)) && (straight = IsCurveStraight( n, tolerance, i, TRUE ));
  478.          i++);
  479.  
  480.     if (straight && n->strV0 && n->strVn)
  481.     n->flatU = TRUE;
  482.  
  483.     straight = TRUE;
  484.     if ( (! n->flatV) && (n->strU0) && (n->strUn) )
  485.     for (i = 1L;
  486.          (i < maxU(n)) && (straight = IsCurveStraight( n, tolerance, i, FALSE ));
  487.          i++);
  488.  
  489.     if (straight && n->strU0 && n->strUn)
  490.     n->flatV = TRUE;
  491.  
  492.     if ( (! n->flatV) || (! n->flatU) )
  493.     return( FALSE );
  494.  
  495.     /* The surface can pass the above tests but still be twisted. */
  496.  
  497.     ScreenProject( &(n->points[0L][0L]),        &cp00 );
  498.     ScreenProject( &(n->points[0L][maxU(n)]),        &cp0n );
  499.     ScreenProject( &(n->points[maxV(n)][0L]),        &cpn0 );
  500.     ScreenProject( &(n->points[maxV(n)][maxU(n)]),  &cpnn );
  501.  
  502.     (void) V3Sub( &cp0n, &cp00, &cp0n );    /* Make edges into vectors */
  503.  
  504.     (void) V3Sub( &cpn0, &cp00, &cpn0 );
  505.  
  506.     /*
  507.      * Compute the plane equation from two adjacent sides, and
  508.      * measure the distance from the far point to the plane.  If it's
  509.      * larger than tolerance, the surface is twisted.
  510.      */
  511.  
  512.     (void) V3Cross( &cpn0, &cp0n, &planeEqn );
  513.  
  514.     (void) V3Normalize( &planeEqn );    /* Normalize to keep adds in sync w/ mults */
  515.  
  516.     d = V3Dot( &planeEqn, &cp00 );
  517.     dist = fabs( V3Dot( &planeEqn, &cpnn ) - d );
  518.  
  519.     if ( dist > tolerance ) /* Surface is twisted */
  520.     return( FALSE );
  521.     else
  522.     return( TRUE );
  523. }
  524.  
  525. /*
  526.  * Turn a sufficiently flat surface into triangles.
  527.  */
  528. static void
  529. EmitTriangles( NurbSurface * n )
  530. {
  531.     Point3 vecnn, vec0n;        /* Diagonal vectors */
  532.     double len2nn, len20n;        /* Diagonal lengths squared */
  533.     double u0, un, v0, vn;        /* Texture coords;
  534.  
  535.     /*
  536.      * Measure the distance along the two diagonals to decide the best
  537.      * way to cut the rectangle into triangles.
  538.      */
  539.  
  540.     (void) V3Sub( &n->c00.point, &n->cnn.point, &vecnn );
  541.     (void) V3Sub( &n->c0n.point, &n->cn0.point, &vec0n );
  542.  
  543.     len2nn = V3SquaredLength( &vecnn ); /* Use these to reject triangles */
  544.     len20n = V3SquaredLength( &vec0n ); /* that are too small to render */
  545.  
  546.     if (MAX(len2nn, len20n) < EPSILON)
  547.     return;                /* Triangles are too small to render */
  548.  
  549.     /*
  550.      * Assign the texture coordinates
  551.      */
  552.     u0 = n->kvU[n->orderU-1L];
  553.     un = n->kvU[n->numU];
  554.     v0 = n->kvV[n->orderV-1L];
  555.     vn = n->kvV[n->numV];
  556.     n->c00.u = u0; n->c00.v = v0;
  557.     n->c0n.u = un; n->c0n.v = v0;
  558.     n->cn0.u = u0; n->cn0.v = vn;
  559.     n->cnn.u = un; n->cnn.v = vn;
  560.  
  561.     /*
  562.      * If any normals are sick, fix them now.
  563.      */
  564.     if ((n->c00.normLen == 0.0) || (n->cnn.normLen == 0.0) || (n->cn0.normLen == 0.0))
  565.     FixNormals( &(n->c00), &(n->cnn), &(n->cn0) );
  566.     if (n->c0n.normLen == 0.0)
  567.     FixNormals( &(n->c00), &(n->c0n), &(n->cnn) );
  568.  
  569.     if ( len2nn < len20n )
  570.     {
  571.     (*DrawTriangle)( &n->c00, &n->cnn, &n->cn0 );
  572.     (*DrawTriangle)( &n->c00, &n->c0n, &n->cnn );
  573.     }
  574.     else
  575.     {
  576.     (*DrawTriangle)( &n->c0n, &n->cnn, &n->cn0 );
  577.     (*DrawTriangle)( &n->c0n, &n->cn0, &n->c00 );
  578.     }
  579. }
  580.  
  581. /*
  582.  * The recursive subdivision algorithm.     Test if the surface is flat.
  583.  * If so, split it into triangles.  Otherwise, split it into two halves,
  584.  * and invoke the procedure on each half.
  585.  */
  586. static void
  587. DoSubdivision( NurbSurface * n, double tolerance, Boolean dirflag, long level )
  588. {
  589.     NurbSurface left, right;    /* ...or top or bottom. Whatever spins your wheels. */
  590.  
  591.     if (TestFlat( n, tolerance ))
  592.     {
  593.     EmitTriangles( n );
  594.     }
  595.     else
  596.     {
  597.     if ( ((! n->flatV) && (! n->flatU)) || ((n->flatV) && (n->flatU)) )
  598.         dirflag = ! dirflag;    /* If twisted or curved in both directions, */
  599.     else                /* then alternate subdivision direction */
  600.     {
  601.         if (n->flatU)        /* Only split in directions that aren't flat */
  602.         dirflag = FALSE;
  603.         else
  604.         dirflag = TRUE;
  605.     }
  606.     SplitSurface( n, &left, &right, dirflag );
  607.     DoSubdivision( &left, tolerance, dirflag, level + 1L );
  608.     DoSubdivision( &right, tolerance, dirflag, level + 1L );
  609.     FreeNurb( &left );
  610.     FreeNurb( &right );        /* Deallocate surfaces made by SplitSurface */
  611.     }
  612. }
  613.  
  614. /*
  615.  * Main entry point for subdivision */
  616. void
  617. DrawSubdivision( NurbSurface * surf )
  618. {
  619.     surf->flatV = FALSE;
  620.     surf->flatU = FALSE;
  621.     surf->strU0 = FALSE;
  622.     surf->strUn = FALSE;
  623.     surf->strV0 = FALSE;
  624.     surf->strVn = FALSE;
  625.  
  626.     /*
  627.      * Initialize the projected corners of the surface
  628.      * and the normals.
  629.      */
  630.     DIVW( &(surf->points[0L][0L]),             &surf->c00.point );
  631.     DIVW( &(surf->points[0L][surf->numU-1L]),         &surf->c0n.point );
  632.     DIVW( &(surf->points[surf->numV-1L][0L]),         &surf->cn0.point );
  633.     DIVW( &(surf->points[surf->numV-1L][surf->numU-1L]), &surf->cnn.point );
  634.  
  635.     GetNormal( surf, 0L, 0L );
  636.     GetNormal( surf, 0L, maxU(surf) );
  637.     GetNormal( surf, maxV(surf), 0L );
  638.     GetNormal( surf, maxV(surf), maxU(surf) );
  639.  
  640.     DoSubdivision( surf, SubdivTolerance, TRUE, 0L );
  641.     /* Note surf is deallocated by the subdivision process */
  642. }
  643.